
## 特征提取和归一化

# ```
# -*- coding: utf-8 -*-
"""
Created on Sun Sep 21 21:47:21 2025

@author: JIAMIN
"""

import numpy as np
import pandas as pd
from scipy.signal import welch, stft
from scipy.stats import kurtosis, skew
import pywt
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
import seaborn as sns


# 中文显示设置
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False
import warnings
warnings.filterwarnings('ignore')

# ==============================================================================
# 数据加载与预览
# ==============================================================================

# 加载源域数据
long_table = pd.read_csv(r'E:\研究生\数学建模\2025研究生题\新\expanded_long_table.csv')
src_long = long_table
print("源域数据预览:")
print(src_long.head())
print("\n源域数据列名:", src_long.columns.tolist())
print("\n源域数据基本信息:")
print(src_long.info())

# ==============================================================================
# 核心特征提取函数
# ==============================================================================

def time_features(x):
    """时域统计特征"""
    x = np.asarray(x).ravel()
    if len(x) == 0: return {}
    rms = np.sqrt(np.mean(x**2))
    mean_abs = np.mean(np.abs(x))
    sqr_mean = np.mean(np.sqrt(np.abs(x)))
    peak = np.max(np.abs(x))
    feats = {
        "mean": np.mean(x), "std": np.std(x), "var": np.var(x), "rms": rms, "peak": peak, "p2p": np.ptp(x),
        "mean_abs": mean_abs, "shape_factor": rms / (mean_abs + 1e-12), "crest_factor": peak / (rms + 1e-12),
        "impulse_factor": peak / (mean_abs + 1e-12), "margin_factor": peak / (sqr_mean**2 + 1e-12),
        "clearance_factor": peak / (sqr_mean + 1e-12), "kurtosis": kurtosis(x), "skewness": skew(x),
    }
    return feats

def freq_features(x, fs, fr):
    """频域特征（转速归一化）"""
    f, Pxx = welch(x, fs=fs, nperseg=2048)
    Pxx = Pxx / np.sum(Pxx)
    f_norm = f / fr
    feats = {
        "spec_centroid": np.sum(f_norm * Pxx),
        "spec_bandwidth": np.sqrt(np.sum(((f_norm - np.mean(f_norm))**2) * Pxx)),
        "spec_skewness": skew(Pxx), "spec_kurtosis": kurtosis(Pxx),
        "spec_entropy": -np.sum(Pxx * np.log(Pxx + 1e-12)),
    }
    bands = [(0.8,1.2),(1.8,2.2),(2.8,3.2),(4.5,5.5)]  #4 个频率区间
    for i,(lo,hi) in enumerate(bands,1):
        mask = (f_norm>=lo)&(f_norm<=hi)
        feats[f"band_energy_{i}"] = np.sum(Pxx[mask])
    return feats

def tf_features(x, fs, fr):
    """时频域特征（STFT + 小波能量）"""
    f,t,Zxx = stft(x, fs=fs, nperseg=1024)
    power = np.abs(Zxx)**2
    power = power / (np.sum(power) + 1e-12)
    f_norm = f / fr
    feats = {
        "tf_entropy": -np.sum(power * np.log(power + 1e-12)),
        "tf_mean_freq": np.sum(np.mean(power,axis=1)*f_norm),
    }
    coeffs = pywt.wavedec(x, 'db4', level=4)
    energy = np.array([np.sum(c**2) for c in coeffs])
    energy_ratio = energy / (np.sum(energy)+1e-12)
    for i,e in enumerate(energy_ratio):
        feats[f"wavelet_energy_{i}"] = e
    return feats

def extract_features_from_long(long_table, signal_col="DE_time", fs=32000):
    """从长表格数据中提取特征（带标签）"""
    feature_rows, labels = [], []
    for fid, group in long_table.groupby("file"):
        x = group[signal_col].dropna().values
        if len(x) == 0: continue
        rpm = group["RPM"].iloc[0]
        fr = (rpm/60.0) if rpm and rpm>0 else 1.0
        feats = {}
        feats.update(time_features(x))
        feats.update(freq_features(x, fs, fr))
        feats.update(tf_features(x, fs, fr))
        feats["file"] = fid
        feature_rows.append(feats)
        labels.append(group["status"].iloc[0])
    X = pd.DataFrame(feature_rows).set_index("file")
    y = pd.Series(labels, index=X.index, name="status")
    return X, y

def extract_features_from_target_long(tgt_long, signal_col="Xtime", fs=32000):
    """从目标域长表格数据中提取特征（不带标签）"""
    feature_rows = []
    for fid, group in tgt_long.groupby("file"):
        x = group[signal_col].dropna().values
        if len(x) == 0: continue
        rpm = group["RPM"].iloc[0]
        fr = (rpm/60.0) if rpm and rpm>0 else 1.0
        feats = {}
        feats.update(time_features(x))
        feats.update(freq_features(x, fs, fr))
        feats.update(tf_features(x, fs, fr))
        feats["file"] = fid
        feature_rows.append(feats)
    X = pd.DataFrame(feature_rows).set_index("file")
    return X

# ==============================================================================
# 执行特征提取
# ==============================================================================

# 源域特征提取
X_src_DE, y_src = extract_features_from_long(long_table, signal_col="DE_time", fs=32000)
X_src_FE, _     = extract_features_from_long(long_table, signal_col="FE_time", fs=32000)
X_src_BA, _     = extract_features_from_long(long_table, signal_col="BA_time", fs=32000)

print("\n源域特征提取结果:")
print(f"特征矩阵形状 (DE): {X_src_DE.shape}")
print(f"标签形状: {y_src.shape}")
print("\n源域标签分布:")
print(y_src.value_counts())


# ==============================================================================
# 目标域处理与可视化
# ==============================================================================

try:
    # 加载并处理目标域数据
    tgt_long = pd.read_csv(r'E:\研究生\数学建模\2025研究生题\新\目标文件\target_long_table.csv')
    print("\n目标域数据加载成功")
    print("目标域数据列名:", tgt_long.columns.tolist())
    
    X_tgt = extract_features_from_target_long(tgt_long, signal_col="Xtime", fs=32000)
    print("\n目标域特征提取结果:")
    print(f"特征矩阵形状: {X_tgt.shape}")

    # 确保特征名对齐
    common_features = sorted(set(X_src_DE.columns) & set(X_tgt.columns))#找源域（X_src_DE）和目标域（X_tgt）共有的特征名
    print(f"共有 {len(common_features)} 个特征可比较")
    print("共有的特征名称:", common_features)
    
    # 特征标准化处理
    all_data = pd.concat([X_src_DE, X_src_FE, X_src_BA, X_tgt], axis=0)
    scaler = StandardScaler()
    scaler.fit(all_data.values)

    #对每个特征矩阵进行标准化转换
    X_src_DE_norm = pd.DataFrame(scaler.transform(X_src_DE), index=X_src_DE.index, columns=X_src_DE.columns)
    X_src_FE_norm = pd.DataFrame(scaler.transform(X_src_FE), index=X_src_FE.index, columns=X_src_FE.columns)
    X_src_BA_norm = pd.DataFrame(scaler.transform(X_src_BA), index=X_src_BA.index, columns=X_src_BA.columns)
    X_tgt_norm    = pd.DataFrame(scaler.transform(X_tgt),    index=X_tgt.index,    columns=X_tgt.columns)


    
    print("\n特征标准化完成")

    import seaborn as sns

    # ==========================================================================
    # 可视化 1: 原始特征比较 (散点图 + 箱线图 并排)
    # ==========================================================================
    print("\n--- 开始生成原始特征比较图 (散点图 & 箱线图) ---")
    for feat in common_features:
        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 5))
        fig.suptitle(f"特征比较: {feat}", fontsize=16)

        scatter_data = [
            X_src_DE[feat].dropna().values,
            X_src_FE[feat].dropna().values,
            X_src_BA[feat].dropna().values,
            X_tgt[feat].dropna().values
        ]
        labels = ["DE_time", "FE_time", "BA_time", "Xtime"]

        # 1. 散点分布图 (stripplot)
        sns.stripplot(
            data=scatter_data,
            ax=ax1,
            palette="Set2",
            jitter=True,  # 增加抖动避免重叠
            size=4
        )
        ax1.set_xticks(range(4))
        ax1.set_xticklabels(labels)
        ax1.set_ylabel(feat)
        ax1.set_title("散点分布图")
        ax1.grid(alpha=0.3)

        # 2. 箱线图
        ax2.boxplot(scatter_data, labels=labels, patch_artist=True,
                    boxprops=dict(facecolor="lightblue", color="blue"),
                    medianprops=dict(color="red", linewidth=2),
                    whiskerprops=dict(color="blue"), capprops=dict(color="blue"))
        ax2.set_ylabel(feat)
        ax2.set_title("箱线图分布")
        ax2.grid(alpha=0.3, linestyle="--")

        plt.tight_layout(rect=[0, 0, 1, 0.95])
        plt.show()

    # ==========================================================================
    # 可视化 2: 归一化特征比较 (散点图替代小提琴图)
    # ==========================================================================
    print("\n--- 开始生成归一化特征比较图 (散点图) ---")
    for feat in common_features:
        plt.figure(figsize=(8, 5))

        scatter_data = [
            X_src_DE_norm[feat].dropna().values,
            X_src_FE_norm[feat].dropna().values,
            X_src_BA_norm[feat].dropna().values,
            X_tgt_norm[feat].dropna().values
        ]
        labels = ["DE_time", "FE_time", "BA_time", "Xtime"]

        # --- 散点图 (stripplot) ---
        sns.stripplot(
            data=scatter_data,
            palette="Set3",
            jitter=True,
            size=4
        )

        plt.xticks(range(4), labels)
        plt.ylabel(feat)
        plt.title(f"特征分布比较 (归一化, Scatter): {feat}")
        plt.grid(alpha=0.3, linestyle="--")
        plt.tight_layout()
        plt.show()

except FileNotFoundError:
    print("\n未找到目标域数据文件 'target_long_table.csv'，跳过目标域特征提取与比较。")
except KeyError as e:
    print(f"\n错误：数据中不存在{e}列，请检查列名是否正确。")
# ```

